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Simulation of inviscid compressible multi-phase 
flow with condensation 


Philip Kelleners f 


Condensation of vapours in rapid expansions of compressible gases is investigated. In 
the case of high temperature gradients the condensation will start at conditions well away 
from thermodynamic equilibrium of the fluid. In those cases homogeneous condensation 
is dominant over heterogeneous condensation. The present work is concerned with de- 
velopment of a simulation tool for computation of high speed compressible flows with 
homogeneous condensation. The resulting flow solver should preferably be accurate and 
robust to be used for simulation of industrial flows in general geometries. 


1. Introduction 

A substance below its critical temperature can be present in either gaseous or liquid 
phase, depending on the pressure, and is referred to as a vapour. Vapours present in 
a mixture of gases and vapours, when subjected to expansion can condensate and form 
liquid droplets. This phenomenon is observed in e.g. aircraft tip vortices and in industrial 
flows like steam turbines. Condensation in flows of gas mixtures at high speed has been 
investigated by, among others; Wegener (1969), Hill (1966), Campbell (1989), Schnerr et 
al. Schnerr (1996), Dohrmann (1989), Mundinger (1994), Adam (1996) and van Dongen 
et al. Luijten (1999), Luijten (1998), Prast (1997) and Lamanna (2000). Expansion in 
nozzles of gases to supersonic speeds has often been used to investigate the physics of 
condensation. Condensation in the flow around airfoil sections and in steam turbines has 
been investigated to a large extent. At the University of Twente a numerical tool has 
been developed to simulate transonic flows with condensation in confined geometries. 
The solver operates on the basis of a finite volume method using unstructured meshes. 
It has been observed that results obtained with the solver are very sensitive to accurate 
shock prediction, and fine shock resolution in the flow field, especially in cases of strong 
interaction between the gasdynamic shock and the condensation process. The focus of the 
present work is to improve the accuracy and robustness of the flow solver by improving 
solid wall boundary treatment and spatial reconstruction for simulations with second 
order spatial accuracy. 


2. Physics of Condensation during Rapid Expansion 

Below its critical temperature, a fluid can be in gaseous or liquid phase. The thermo- 
dynamical region of coexistence of vapour and liquid in equilibrium in bulk substances 
is given by the Clausius-Clapeyron relation. For rapid expansions of vapour this coexis- 
tence region is passed without the fluid attaining equilibrium. The vapour is saturated 
expressed by the saturation ratio; 
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Figure 1 . Expansion of air-water mixture in nozzle-S2, partial vapour pressures and 

saturation ratio s (dashed) 


where n v and n v , eq are the number of moles of the vapour in the actual mixture and 
the mixture in equilibrium, respectively. For pressures as high as atmospheric pressure 
equation 2.1 can be replaced by the more conventional expression: 


where p v is the actual partial vapour pressure in the mixture and p v ,eq is the equilibrium 
partial vapour pressure at the same thermodynamic conditions. s,p v and p v _ eq are plotted 
in figure 1 for the expansion of an air- water mixture, as functions of temperature. 

The curve for the coexistence of liquid and vapour (p Vt eq) divides the plane given by 
pressure and temperature into two regions. For pressures higher than the coexistence 
pressure the substance under consideration (water) is present in liquid form while in 
equilibrium state. For lower pressures the substance will be present as water-vapour 
while in equilibrium. The expansion in the nozzle, as given by the partial vapour pressure 
p v , is so quick that the water- vapour will not immediately condense under equilibrium 
conditions. This is the case as the characteristic time of the gasdynamic flow is much 
smaller than the time needed to form the first onsets of the new liquid phase. For the 
nozzle flow at hand this is due to the high-cooling rate in the nozzle. So the water-vapour 
expands further, driving the air-water-vapour mixture well away from equilibrium, as 
indicated by the saturation ratio which attains values as high as 40. In case s > 1 the fluid 
is said to be super-saturated. Formation of small liquid clusters, nuclei, at high super- 
saturation is the first stage of the condensation process that starts in order to reestablish 
equilibrium. On these newly formed nuclei, the super-saturated vapour condenses as a 
second step until equilibrium is reached, the process of droplet growth. This can be seen 
in figure 1 by the decrease in partial vapour pressure and the increase in temperature. 
With so much liquid already present, the remainder of the condensation process due 
to droplet growth takes place very close to thermodynamic equilibrium of the water, 
indicated by the saturation close to one, for the remainder of the expansion. 
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The process described above is the process of homogeneous condensation. In contrast, 
the process of heterogeneous condensation is droplet growth on foreign, already present 
particles. However, for high cooling speeds and consequentially high super-saturation 
the number of new nuclei formed exceeds any realistic number of foreign particles by 
several orders of magnitude. For the present applications only homogeneous condensation 
processes are of importance. 

An effect resulting from condensation is the release of latent heat during the con- 
densation process, indicated by the increase in temperature during condensation in the 
previous example of nozzle flow. The latent heat L is defined as 

L = h v — hi (2.3) 

with h v and hi the enthalpy of the gaseous and liquid phase, respectively. The latent 
heat is the enthalpy needed to evaporate a unit mass of liquid. It is a material property. 


3. Nucleation and Droplet Growth models 

From the previous treatment of condensation during rapid expansion, it is observed 
that the condensation process consists of two consecutive stages. The first one is the 
formation of liquid nuclei, nucleation, the second one is the condensation of vapour 
molecules on the already present nuclei, making these grow in the process of droplet 
growth. In the following, a brief treatment is given of the physical background of both 
nucleation and droplet growth together with the presentation of the models used to 
simulate these processes. 

3.1. Nucleation 

Under super-saturation vapour molecules can condense on liquid already present or on 
foreign objects. In the absence of both, the vapour molecules can form small clusters. As 
a result of the clustering, an additional phase needs to be formed, the interface between 
the liquid inside the cluster and the gas outside of the cluster. The interface can be re- 
garded as infinitely thin. At the interface there is surface tension. Thus the creation of an 
interface requires energy. For very small clusters, the surface effects are dominant over 
volume related effects. As a result the formation of the interface represents an energy 
barrier in the formation process of the nucleus. If the energy involved in the clustering 
of the vapour molecules is less than required in the formation of the interface surface, 
the cluster will disintegrate immediately following its formation. Therefore at near equi- 
librium conditions, super-saturations close to one, it is highly unlikely, that a realistic 
number of stable nuclei in a volume at macro scale will be formed although clusters 
are constantly formed and falling apart. The previous notion results in formulation of 
formation enthalpy of critically sized stable nuclei, and the number densities in which 
these stable nuclei are likely to come into existence. The creation enthalpy of one liquid 
droplet under equilibrium conditions is given by Dohrmann (1989): 

AG = 4nr 2 a — nkT ln(s) (3-1) 

where, AG is the change in Gibbs enthalpy, r is the radius of the droplet, a is the droplet 
surface tension in a plane with no curvature, n is the number of molecules in the droplet, 
k is the Boltzmann constant, and T is the temperature. In accordance with classical 
nucleation theory, we note that a stable nucleus will be created when the function for 
AG attains a local maximum. AG at this maximum is: 
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AG = ^irr* 2 a (3.2) 

with r* , the radius of this critical sized droplet being given by: 

= piR v T\n{s) (3 ' 3) 

where pi is the density of the liquid in the droplet, and R v is the gas constant of the 
vapour. The number of droplets created is given by the nucleation rate. For the nucleation 
rate several models are available, Classical Nucleation Theory CNT-model (Wegener 
1969) and the Internally Consistent Classical Theory, ICCT-model (Luijten 1998). In 
the present work, the CNT-model is applied: 


pf- 1 ? , ^ 2 ) P.4) 

Pi V 7rm 3 Y 3 mpf R 3 T 3 In (s) ) 

where p v is the density of the vapour and m is the mass of one vapour molecule. 

3.2. Droplet Growth 

The droplet growth model to be applied depends on the regime in which droplet growth 
takes place. For pressures of the order of atmospheric pressure droplet growth is based 
on a balance between condensation of vapour molecules onto droplets, and evaporation 
of vapour molecules from the droplet. For pressures 1 to 2 orders of magnitude higher, 
droplet growth is diffusion-controlled. In the present work, high pressure effects are not 
taken into account. Therefore the Hertz-Knudsen droplet growth model (Hill 1966) can 
be applied. The droplet growth rate is: 


dr _ a ( p v Ps, r__\ ^ ^ 

dt pi \vZ27r TRjT \J2 i iR v T d j 

where a is an accommodation coefficient usually taken equal to one, T<j is the droplet 
temperature in the present work equal to the surrounding gas temperature and is 
the partial super-saturated vapour pressure over a curved radius of curvature r given by: 


Ps,r = Pv,eq exp 


2a \ 

piR v Tr h i J 


(3.6) 


where rhi is the Hill droplet radius. This radius is computed as an averaged radius over 
the complete population of droplets at hand. This radius, of course, depends on the 
distribution function of the droplets. 


4. Description of the Liquid Phase 

Until now only, the creation and growth rate of a single droplet, and the creation rate, 
nucleation rate, of new droplets has been treated. Next to that, the size and form of 
individual droplets and the size and distribution of the droplet population need to be 
computed. As the number of droplets in the flow cases at hand typically range between 
10 12 and 10 22 number of droplets per unit mass, it is impossible to describe individual 
droplets. One alternative approach is to define classes of droplets of similar radius, the 
class-model. This approach naturally extends to a description of the distribution of the 
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Figure 2. Liquid droplets in control volume 


droplets. For the nozzle flow problems at hand, the typical range of radii of the droplets 
will be between 2. ICC 10 and 1.10“ This would require an impractical high number of 
droplet-classes, making the class-model computationally very expensive. Hill’s method 
of moments is computationally less intensive at the cost of some of the resolution of the 
droplet distribution. As for the problems at hand, the total amount of liquid generated 
and the strong interaction between the gasdynamic flow and the condensation process is 
of primary interest, the loss of resolution of the droplet distribution is no severe penalty. 

4.1. Hill’s Method of Moments 

Consider a control volume as depicted in figure 2. In the control volume a mixture of an 
ambient gas, and a vapour in its gaseous and its liquid state is present. A first assumption 
is that all liquid is present in the form of spherical droplets. Conservation of the mass in 
the control volume gives: 


M = M a + M v + Mi (4.1) 

where M is the total mass in the control volume, M a is the mass of the ambient gas, 
M v is the mass of the condensible substance in gaseous form and Mi is the mass of the 
same substance in liquid form. The ambient air is assumed to be permanently above its 
critical temperature. The dimensionless liquid mass fraction g is defined as: 


(4.2) 


_ Mi _ Pl Vi _ Vi 
9 M P V Pl pV 
With the final notation emphasis is placed on the next assumption; the liquid mass 
present in the control volume is regarded to be incompressible. Therefore pi is taken 
constant. The liquid mass fraction is a function of time. The density of the complete 
mixture; gas, vapour and liquid p, is dependent of time, as are the volume occupied by 
the liquid V) and the total control volume V. This gives: 

^ _ Vi{t) 

9 pl P (t)v(ty 

The volume occupied by the liquid V) is given by the following integral: 


t 

Vi(t) = J ^7rr 3 (f,T)J(r)C(r)dr (4.3) 

0 

where J(r) is the nucleation rate per unit volume and V{r) is the size of the control 
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volume, at the moment of nucleation. The volume occupied by the droplets is given 
by |7rr 3 (t,r). The function r(t,r) describes the radius of the droplet, from its time of 
creation, r(ro, To) where this radius is equal to the critical radius as given in equation 3.3, 
r(r 0 ,r 0 ) = r* , up to the present radius of the droplet at time t. Under the integral the 
product J(tq)V{tq) acts as a weight function for the contribution of the liquid volume 
of the droplets created at time tq to the total amount of liquid formed at time t. The 
liquid mass fraction now becomes: 


9{t) = Pi 


t 


I 


|7r r 3 (f, t)J(t)V (r)dr 


p{t)V{t) 


If the considered lump of mass, liquid and gas-vapour mixture, is at rest, or flows along 
a stream-line, the fraction is fixed. This is identical to the assumption that the 

droplets do not move relative to the surrounding gas mixture. This is known as the 
no-slip condition in Hill’s Method of Moments. This can be expressed as: 


dM = 0 => M = p(t)V(t) = p(t)V(t) 

The liquid mass fraction can be rewritten (Hagmeijer 2001): 

g(t)=pi j ^7rr 3 (t,r)^ydr. 
o 

However, the function r(f, r) is not readily available in closed form. The droplet growth 
rate % is known in closed form. Therefore g(t) is differentiated with respect to time: 


dg d 
dt ^ 1 dt 


t. 


o 



J(t) 

p(t) 


dr 


= Pi 



dr(t , r) J(r) 
dt p(r) 


dr + pi-7rr 3 (t,t) 


M 

p(t) 


4 

= Pi -7-c + Zpi 


P(t) 



dr(t,r ) J(t) 

; 

dt p{r) 


In this analysis, it is assumed, that at the initial time t = 0 there is no liquid present 
in the control volume. A third assumption is that the present droplet growth rate is 
independent of the present radius of the droplet. This assumption is certainly not true 
for very small droplets because of the dependency of the surface tension on the droplet 
radius, however the assumption is valid for droplets with larger radii. 


dr , dr , , 

H* dt (r(t - T)) (4 - 4) 

This implies that, the present value of the droplet growth rate, is independent of its 
history. So it can be written: 


dr dr 
dt dt 


(4.5) 





Inspection of the last moment, Q o, shows that there is no longer a functional dependence 
of t in the integral. So the differentiation of this last liquid moment can be expected to 
be different from the higher liquid moments: 


dQ o = d_ ( f J(t) \ = J(t) 
dt dt yj p(r) T J p(t) ' 
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This last result is fortunate, as it shows, that the logical expansion defined by the recur- 
rence relation for the liquid moments, is not continued for n < 0. So there is no closure 
problem for the solution of the set of moments. The system of ordinary differential equa- 
tions, to determine the liquid moment Q 3 now becomes: 


dQs 

dt 

dQ2 

dt 

dQ\ 

dt 

dQp 

dt 


" Pit) 

Pit) 

r * 

Pit) 

m 

pit) 


- 2 S< w > 


In this system the variables represent the following quantities: 

• Q 3 is proportional to the sum of all droplet volumes, 

• Q 2 is proportional to the sum of all droplet surface areas, 

• Qi is proportional to the sum of all droplet radii, 

• Q 0 is the present number of droplets. 

Now it is assumed that the droplet growth rate %{t) is given only by the present local 
thermodynamic and chemical state, and not by the spatial gradients of r in the flow 
domain. In this case, = 0. By addition of the product of vector Q and the continuity 
equation for mass, the system of ordinary differential equations above can be rewritten 
in strong-conservation form: 


dpQ 3 

dt 

d P Q 2 

dt 

dpQi 

dt 

dpQo 

dt 


dpQ 3 Uj 

dxi 

dpQ 2 Uj 

dxi 

dpQiUi 

dxi 

dpQpUj 

dxi 


= r* 3 J(t) + 3 — (t)pQ 2 
= r* 2 J(t) + 2^(t)pQ 1 
= r*J(t)+^(t)pQ 0 

at 

= m 


4 ipQsit)) 

git ) = pi~ n 

3 p 


Where ay and iq are the spatial coordinates, and the velocity components in spatial 
direction respectively, r * , J and % (t) are given by the laws formulated for the creation 
and growth of a single droplet. The above set of equations gives the integral properties 
of the droplet distribution. So details of individual droplets are not available. To be able 
to compute a radius needed for droplet growth laws, Hill defined (Hill 1966) an averaged 
radius as: 



Thl = 


(4.6) 
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Reason for application of Q 2 in calculation of ?>,; , is the dominance of surface effects in 
droplet growth, hence the application of the integral surface quantity Q 2 - 


5. Thermodynamic State of Gas Liquid Mixtures 

The change of phase from vapour to liquid releases latent heat to the surrounding 
mixture. Due to the condensation, the fractions of the gases in the mixtures change. 
Both processes result in a change of the thermodynamic state of the mixture. To model 
these changes regard the enthalpy of the mixture: 


M h — i\4 a h a M v h v Mi hi 

where h = H /M is the specific enthalpy. By application of the previously defined dimen- 
sionless mass fraction g = Ff : 

h = (1 gmax)ha ( gmax g^)h v ghi . 

Application of the definition of the latent heat L = h v — hi and the definition c p = ( ) p g 
results in: 

dL 

Cp = (1 gmax)c Pa + ( gmax ~ g) c p v A 17(Cp„ ~ ffj' . (5-1) 

Under the assumption — — ss — , a similar expression can be derived for the isochoric 
specific heat coefficient c v = (§^) Vg : 

dL 

Cy = (1 ~ gmax)C Va + ( gmax g) ( 'V v H - ff( c p„ — )• (5-2) 

With these relations for the specific heat coefficients, the gas constant R m ix and the 
dimensionless ratio 7 m ix can be computed similar to the case of an ideal gas, R m ix (5, L) = 
c p — c v and 7 m i X (giL) = 'ff . Both R and 7 are now functions depending on the mass 
fraction g and the relation for the latent heat L. Derived quantities like pressure p and 
speed of sound c can be shown (Mundinger 1994) to be: 

p=pR(g,i)T c 2 = 7 ( 5 , l )-. 

p 


6. Model of Inviscid Flow with Condensation 

Based on Reynolds numbers for typical flow problem of interest the flow is assumed 
to be inviscid. The Euler equations are used to describe the flow. Together with Hill’s 
Method of Moments, the models for nucleation and droplet growth and an equation of 
state they form a closed set of equations: 

| -J<p'dv+ j £{<£')■ ds = J W(<p')dv ( 6 . 1 ) 

V S=3V V 

where <j>' is the vector with the conserved quantities, F(<p') the flux vector, and W{4>') is 
the source term: 
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puu + Ip 
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puH 


0 

pQz 

m') = 

puQ 3 

W{$) = 

r* 3 J(t)+3%(t) P Q2 

pQ-2 


puQ 2 


r* 2 J(t) + 2%(t) P Q 1 

pQi 


puQ i 


r*J(t) + %(t) P Qo 

_ pQo _ 


puQo 


l At) \ 


where p , u,p are the density, the velocity and the pressure of the mixture. E = e+ |||u|| 2 
and H = E + 2 are the total energy and the total enthalpy of the mixture. Qi is the i-th 
moment in Hifi’s Method of Moments. The thermal equation of state reads: 

p = P R(g,L)T (6.3) 

and the caloric equation of state: 

e = c v (g,L)T. (6.4) 


7. Finite Volume Discretization 

With the definition of the control volume averaged (ftp. 



Vi 


equation 6. lis semi-discretized as: 

^ Vi + • s itj = w(4>)Vi 

3 = 1 
n 

The summation of discrete fluxes fiAj) ' &i,j is computed using an edge-based data 

3 - 1 

structure as described in Jameson (1986), Barth (1989) and Barth (1994). For every 
edge the indices of the 2 control volumes connected by the edge are stored, as well as 
the three spatial components of the surface normal vector of the surface between the two 
control volumes. The length of the surface normal vector is equal to the magnitude of the 
surface area. This edge-based approach allows for a cell-centered or a vertex-centered 
use of the original mesh, without significant changes to the flux computation algorithm. 
The fluxes are computed, by solving a local one-dimensional Riemann problem for every 
surface using an approximate Riemann solver. This is implemented in the flow solver for 
edges oriented in every possible direction in three-dimensional space. Two-dimensional 
problems with meshes in two independent directions can be regarded as a subclass of the 
full three-dimensional problem. The same holds for quasi one-dimensional or truly one- 
dimensional flow problems. With proper definition of additional boundary surfaces in the 
case of quasi one-dimensional flows, all flows can be computed using the same flow-solver 
regardless of the mesh being fully three-dimensional, two-dimensional or quasi or truly 
one-dimensional. The preprocessor, to the flow solver, generating the edge-based data 
structure from the meshes needs only to insert zeros for the edge-related surface normals 
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in the directions not relevant for the particular mesh. To date, pre-processors have been 
written for meshes consisting of line elements for one-dimensional meshes, triangular ele- 
ments for unstructured two-dimensional meshes, quadrilaterals for structured monoblock 
two-dimensional meshes, and tetrahedral elements for three-dimensional meshes. How- 
ever, the edge-based data structure allows simple extension to multiblock meshes, meshes 
constructed by hexahedra and hybrid meshes. This requires only the pre-processor for 
the particular element(s) to be developed, and minor additions to the domain boundary 
treatment in the flow solver. The flow solver has successfully completed simulations for 
one-, two- and three-dimensional flows. The great advantage of this approach is that 
almost all testing during development can be done for flow problems in one and two- 
dimensions decreasing development time, as formulation of one and two-dimensional flow 
problems and generation of one- and two-dimensional meshes is much less time consum- 
ing than in the full three-dimensional case. There is a penalty in the case of computation 
of one- and two-dimensional flows. For the one-dimensional case the momentum fluxes 
in the two spatial directions perpendicular to the main spatial direction are computed 
but not used. In the case of the Euler-equations this would result in 66% increased com- 
putational expense. In the two-dimensional case this increased computational expense 
is 25%. In the case of the computation of the Euler-equations and Hill’s momentum 
equations, this computational overhead drops to 28% in the one-dimensional case, and 
12% in the two dimensional case. However, the number of grid points for flow prob- 
lems in two-dimensions and one-dimensions is one and two orders lower respectively. So 
the computational overhead is completely negligible next to the enormous advantage of 
simple and faster testing in lower dimensions. 

Eigenvalue-analysis Kelleners (2001) of the complete system of Euler-equations with 
the Hill Momentum equations has shown that the additional eigenvalues due to the Hill 
equations are all real and have value u. So the liquid moments Qo-.Qz are convected 
downstream along streamlines, and no new acoustic waves are introduced. The computa- 
tion of the flux is identical to the case of flow without condensation, with only additional 
transport equations for the liquid moments. 

The presence of the source terms in the Hill momentum equations results in the com- 
plete system of partial differential equations being a stiff system. Wishing to use an ex- 
plicit time-stepping method as a relaxation method to compute a steady state, the fact 
that the system is stiff would lead to an unacceptable small time-step requirement due 
to the source terms. To circumvent this time-step restriction a fractional time-stepping 
method is used as proposed by Oran and Boris, Oran (1987) the differential equation of 
the form: 


-^+m = w{4>) 

the solution procedure is split into: 

(7.1) 

(7.2) 


df* 

dti 

d(j> n+1 

dU 


= -f{<n 

= w(<j>*) 


where the time-step operator for the first step can be any conventional explicit operator, 
e.g. Euler-forward or Runge-Kutta. The time step-operator for the second time step 
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needs to be able to cope with the possible large magnitude of the source w(<f>*), (see 
Oran 1987; Mundinger 1994) 

The two following subsections highlight two topics which require detailed attention in 
the case of development of a spatial second-order accurate finite-volume method. For 
other topics like, flux splitting or flux limiting, see Jameson (1995) and Liou (1996). 

7.1. Higher Order Spatial Reconstruction 

To achieve good accuracy on meshes with moderate vertex densities the computational 
method must be at least second order. Whereas computation of first-order spatial ac- 
curacy allows for the storage location of the cell-averaged value to be located anywhere 
inside the control volume, this point should preferably be located at the geometric center 
of gravity for higher-order spatial accuracy. Use of the center of gravity circumvents se- 
vere penalties to the applicable CFL-number of individual control volumes. Higher-order 
reconstructions of cell-averaged data should conserve the cell-averaged value. In case of 
second-order reconstruction it can be proven that use of the center of gravity satisfies 
this requirement: 


u{x) = u cg + — (a? - x cg ) (7.3) 

where u cg is the cell-averaged value positioned at the center of gravity. u(x) integrated 
over the control volume gives: 

Y J v{x)dv = ^ J U cg dv + ^ J ^ {X - x cg )dv (7.4) 

V V v ~ ° 9 

The first term on the right-hand side gives the cell-averaged value. The second term on 

dU 

the right hand side is zero. The gradient can be moved in front of the integral. The 

-eg 

remaining integral is the static moment about the center of gravity, and necessarily zero. 
For vertex-centered control volumes, the movement of the cell storage location, from 
cell-vertex to center of gravity of the control volume, is largest for control volumes at the 
physical boundaries of the domain. 

7.2. Linear Reconstruction at Domain Boundaries 
The assumption of linear reconstruction in the control volume in the case of second 
order spatial accuracy, requires special attention at the domain boundaries. At solid wall 
boundaries the pressure part of the flux needs to be calculated. Question is, how the 
pressure must be integrated along the solid wall boundary surface for the resulting flux 
integral along the complete boundary of the control volume to be of second-order spatial 
accuracy. It is assumed that fluxes vary linearly over the control volume. Reconstruction 
of the flux at the solid wall boundary is therefore similar to linear reconstruction from 
the cell averaged state of the conservative variables throughout the control volume. The 
latter requires the computation of gradient the ^ at the center of the control volume. 
Green-Gauss reconstruction is often used for computation of this gradient Barth (1994); 
Aftosmis (1995). Green-Gauss reconstruction is easily implemented using the edge-based 
data structure. Starting point is Gauss’ divergence theorem: 


V (4>)dz 


(f>ds 


v 


(7.5) 
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Figure 3. Median dual mesh control volume sharing sides with a physical boundary 
which is rewritten in discrete form: 

= (™) 

i 

As an example the linear distribution of <f> along a domain boundary for a control volume 
being the median dual to a two-dimensional triangular element is given, see figure 3. 
The volume V is the surface area associated with vertex 0 in the case of the median dual 
mesh is: 

V=l\\QXb\\ 

To evaluate equation 7.6, the outward pointing surface vectors and the mean 

values of the quantity need to be determined. 

1 

s_i = -a x e 2 

A , 1 , 

Sj 2 = (gO- gQj X e_ z 

A, l , 

S 3 = A- xSj X e z 
“63 

s _4 = x e_ 2 

where e z = e x x e y . <j>i is computed similar as would be the case for a median dual mesh 
volume embedded entirely in the physical domain. The value for </>j at surfaces s 2 and s 3 
is interpolated: 


< / > s 2 — 2 ^° + ^ 

< j > s 3 = 2(^0 + &) 

At the surfaces and s 4 it is assumed that there exists a linear distribution of f> along 
the surface of the following form: 

<j) = axj ) 0 + (1 - u>)4>i 
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The allowable value of u> is between zero and one. Substitution of s, ; and <j>i in equation 
7.6 ultimately gives: 


VO) = = 

i 

((3a; — |)0o +(|-3w)0i -<j> 2 ) ^ 

((—3a; +§)0o +4>i +( _ |+3o;)02) 

To compute the value of u>, the gradient V (<j>) is also computed assuming a linear distri- 
bution of <f> in the triangle given by the vectors a and b. 



(t>(x) -<j> 0 = a(x - x 0 ) + /3{y - y 0 ) 
The gradient for this linear distribution is: 


VO) = I = 


dy 


Using the data in vertices 1 and 2, <pi,<j )2 a system can be derived, which can be solved 

f Qt \ 

for I o ) ■ The solution of this system is: 


( a \ _ f a x a y \ ( <Pi ~ <Po \ 

\P ) \b x b y J { O - <h ) 


Both expressions for V (< t > ) are equal to one— another. Formulating this identity, and noting 
that it must be true for any value of (po. </>i jOj it can be shown that: 


L) = 


5 

6 


Returning to the assumption that the wall pressure flux is linearly distributed along the 
domain boundary surfaces, in a manner similar as any reconstructed scalar </>, the weights 
| and g can be applied to compute the integral of the pressure over the domain boundary 
surfaces. E.g. for surface s 4 this becomes: 


Ps 4 s 4 = (gPo+ gP2)s 4 . 


It should be stressed that the value of u> depends on the type of element at the domain 
boundary, (e.g. triangle, tetrahedron or brick) and the manner in which the control 
volume is defined (e.g. cell-centered or vertex-centered median dual mesh). 


8. Results 

The simulation tool described above has been used to compute solutions to flow prob- 
lems of internal and external flows, adiabatic flows and flows with condensation. The 
cases presented below have been computed using the following layout of the solution 
algorithm. Node-centered finite-volume scheme, using median dual mesh cells. Explicit 
time stepping of the homogeneous equations as described in equation 7.2 using the second 
order Heun-method. For second order spatial resolution, application of a least-squares 
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algorithm for computation of spatial gradients in the flow field. Reconstruction of limited 
second-order states at the control volume interfaces using minimum-maximum restric- 
tions formulated by Barth (1989), and a limiter function by Venkatakrishnan as presented 
in Aftosmis (1995). Fluxes computed using flux splitting according to Eberle. Compu- 
tation of the time-updates due to the source terms as mentioned in equation 7.2 using 
a fractional time stepping method as presented in Mundinger (1994). Inflow boundaries 
conserving entropy and total enthalpy at infinity upstream for a general equation of state, 
outflow boundaries maintaining downstream static pressure in case of subsonic outflow, 
or simple first-order extrapolation of flow quantities in case of supersonic flow. At solid 
walls application of linear pressure distribution along the outward facing surfaces, as 
presented in section 7.2. Heat effects due to condensation, are taken into account by in- 
clusion of both the gaseous and liquid phase into the energy conservation equation. This 
results in strong two-way coupling of the condensation process and the gasdynamics. 
All computations related to the equation of state, or material properties have been pro- 
grammed into separate library routines. These libraries are called upon by the routines 
solving the conservation equations. This produces additional computational overhead, 
but should allow for quick implementation of different equations of state. In the present 
form of the flow solver only an equation of state for an ideal gas is implemented. 

8.1. Condensation in Nozzle Ba- 120 

Flows with condensation in nozzle Ba-120, designed by Bartlma, were investigated both 
experimentally and numerically by Schnerr and co-workers (see Mundinger 1994). Flows 
for two different humidities are presented in figure 4. The stagnation conditions for both 
flows in figure 4 are: 

s 0 [%] T o [K] Po [pa] 

42.1 298.8 100900 
49.3 297.1 100900 

The expansion of the air-water-vapour mixture gives rise to condensation in the flow 
field downstream of the nozzle throat. This is seen by the steep rise in the nucleation 
rate. Following the creation of the nuclei is the process of droplet growth indicated by 
the rise of the liquid mass fraction 9 . The release of latent heat to the flow, induces 
a shock in both cases. However the shock-strength in the case of the higher humidity is 
much larger, and as a result the nucleation pulse is terminated abruptly. In both cases, 
the continuing expansion due to nozzle divergence results in condensation of almost all of 
the water-vapour at nozzle station x = .15[m]. In the Mach-isoline plot of the flow with 
higher humidity, the reflection of the curved shock-wave from the nozzle wall is nicely 
visible. 


8.2. Condensation in Vortex Flow 

To study the process of condensation in vortical flow, a very slender delta wing has been 
placed inside a tube, see figure 5. At the inlet of the tube a mixture of air-water-vapour 
flows in at Mach 1.55. The humidity of the mixture at stagnation condition is as low as 
1.1%. The wing induces a vortex in the flow field. In the vortex core there is considerable 
loss of total pressure, clearly visible in figure 6. The vortex spirals downstream in the 
tube in a helical form. In the low pressure region in the vortex on the upper side of the 
wing, super-saturation occurs, resulting in large nucleation rate, see figure 7. Note that 
the region with the highest nucleation rates remains restricted to the front half of the 
wing close to the apex. This is a consequence of the droplet growth lowering the levels of 
super-saturation considerably once the first liquid droplets are created in great numbers. 
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Mach, AM=.02 




Figure 4. Flow with condensation in nozzle Bal20, flow from left to right. Humidity at stag- 
nation condition for left figure: so = 42.1%, right figure so = 49.3%. Mach iso-lines in upper 
figures, only iso-lines for M > 1 are drawn. Nucleation rate J, dimensionless pressure ^ , and 
liquid mass fraction ■' — in lower figures. Dashed lines indicate center line, and position of 

nozzle throat 

The dimensionless liquid mass fraction 3 is plotted in figure 8. On the second half of 
the upper side of the wing, almost all of the water is present in liquid form. Beyond the 
trailing edge of the wing, some evaporation of the liquid water occurs, but the main part 
of the water is convected downstream in liquid form. As the vortex is diffusing over the 
cross-sections of the tube, so is the cloud of water droplets, as can be noted from the 
similar shapes of the iso-lines of both the total pressure in figure 6 and the liquid mass 
fraction in figure 8 


9. Conclusions 

A brief treatment has been given on condensation in high speed flows. Because of large 
temperature gradients, the condensation can start well away from thermodynamic equi- 
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Figure 5. Layout of the delta wing, which 
produces the vortex. For clarity only one 
half of the tube wall is shown. Flow from 
upper right corner to lower left corner. 



Figure 6. Total pressure, pt 




Figure 7. Nucleation rate, J plotted on a 
log-scale. 


Figure 8. Liquid mass fraction, g. 




librium, leading to fog-like clouds consisting of many small droplets. Basic nucleation 
and droplet-growth models have been presented. A thorough derivation of the integral 
description of the liquid phase, Hill’s Method of Moments has been given. The com- 
plete model describing inviscid compressible multi-phase flow with condensation and its 
finite-volume discretization have been given. Two aspects of the implementation of the 
finite-volume method for achieving second-order spatial accuracy have been highlighted. 
Results presented for a two-dimensional flow problem, and a fully three-dimensional flow 
problem have been computed using one and the same flow solver. The benefit of develop- 
ing only one single flow solver is clearly felt. From the case of higher humidity in nozzle 
Ba-120, a strong interaction between gasdynamics and condensation can already be seen. 
Higher humidities will lead to unsteady flow. Both experimental and numerical data is 
already available from Schnerr and co-workers. A future extension of the flow solver to 
take into account unsteady flow should be made. For flows with condensation in vortices 
simulations should be performed using the full Navier-Stokes equations, to validate the 
assumption of inviscid flow, made in this work. Improvements must be made with respect 
to accuracy and robustness of the solution algorithm, to make results of the simulation 
tool less dependent on grid quality, but the last problem is universal to CFD-methods. 
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